Prac 5 - IT de Villiers (17502292)¶
In [420]:
#All the necessary imports
%matplotlib inline
import pylab as pl
import numpy as np
from scipy import signal
import IPython.display
pl.rcParams['figure.figsize'] = (9,2)
def setup_plot(title, y_label, x_label):
pl.box(False)
pl.margins(*(pl.array(pl.margins())+0.05))
pl.title(title)
pl.ylabel(y_label)
pl.xlabel(x_label)
In [421]:
def zplane(zeros, poles):
"""
Plot the complex z-plane given zeros and poles.
"""
zeros=np.array(zeros);
poles=np.array(poles);
ax = pl.gca()
# Add unit circle and zero axes
unit_circle = pl.matplotlib.patches.Circle((0,0), radius=1, fill=False,
color='black', ls='solid', alpha=0.6)
ax.add_patch(unit_circle)
pl.axvline(0, color='0.7')
pl.axhline(0, color='0.7')
#Rescale to a nice size
rscale = 1.2 * np.amax(np.concatenate((abs(zeros), abs(poles), [1])))
pl.axis('scaled')
pl.axis([-rscale, rscale, -rscale, rscale])
# Plot the poles and zeros
polesplot = pl.plot(poles.real, poles.imag, 'x', markersize=9)
zerosplot = pl.plot(zeros.real, zeros.imag, 'o', markersize=9, color='none',
markeredgecolor=polesplot[0].get_color(),
)
#Draw overlap text
overlap_txt = []
def draw_overlap_text():
for txt in overlap_txt:
try: txt.remove()
except: txt.set_visible(False)
del overlap_txt[:]
poles_pixel_positions = ax.transData.transform(np.vstack(polesplot[0].get_data()).T)
zeros_pixel_positions = ax.transData.transform(np.vstack(zerosplot[0].get_data()).T)
for (zps_pixels, zps) in [(poles_pixel_positions, poles), (zeros_pixel_positions, zeros)]:
superscript = np.ones(len(zps))
for i in range(len(zps)):
for j in range(i+1,len(zps)):
if superscript[i]!=-1:
if np.all(np.abs(zps_pixels[i] - zps_pixels[j]) < 0.9):
superscript[i]+=1;
superscript[j]=-1;
for i in range(len(zps)):
if superscript[i] > 1:
txt = pl.text(zps[i].real, zps[i].imag,
r'${}^{%d}$'%superscript[i], fontsize=20
)
overlap_txt.append(txt)
draw_overlap_text()
#Reset when zooming
def on_zoom_change(axes): draw_overlap_text()
ax.callbacks.connect('xlim_changed', on_zoom_change)
ax.callbacks.connect('ylim_changed', on_zoom_change)
In [422]:
import os
import urllib
import scipy.io
from scipy.io import wavfile
#Download yesterday.wav from courses.ee.sun.ac.za and return it as a numpy array
def yesterday_wav():
url = 'http://courses.ee.sun.ac.za/Stelsels_en_Seine_414/content/yesterday.wav'
filename = os.path.split(url)[-1]
#Download if path does not already exist
if not os.path.isfile(filename):
urllib.request.urlretrieve(url, filename)
sample_frequency, signal_array = wavfile.read(filename)
#Normalise signal and return
signal_array = signal_array/np.max([np.max(signal_array), -np.min(signal_array)])
return sample_frequency, signal_array
In [423]:
#Example of subplotting
pl.figure(figsize=(11,2))
pl.subplot(1, 2, 1)
setup_plot('plot 1', 'y-label', 'x-label')
pl.stem([1,2,3,4])
pl.subplot(1, 2, 2)
setup_plot('plot 2', 'y-label', 'x-label')
pl.plot([4,3,2,1]);
Question 1: All Zeros Comb filter¶
In [424]:
def sinewave(F,t):
return np.sin(2*np.pi*F*t)
fs, x1 = yesterday_wav()
In [425]:
f0 = 2205
n = np.arange(len(x1))
T = 1/fs
t = n*T
x0 = sinewave(f,t)
x2 = np.abs(x0)
x = x2+x1
In [426]:
X1 = np.abs(np.fft.fft(x1))
X2 = np.abs(np.fft.fft(x2))
X = np.abs(np.fft.fft(x))
k = np.linspace(0, fs/(2*np.pi), X.size)
k2 = np.linspace(0, 1, X.size)
time_range = np.logical_and(t>=3,t<3.005) #in order to display only between 3 and 3.005
In [427]:
pl.figure(figsize=(11,2))
setup_plot('x1(t) - Yesterday - The Beatles', '|x1(t)|', 'time (s)')
pl.plot(t[time_range], x1[time_range]) #snippet between 3s and 3.005s
pl.figure(figsize=(11,2))
setup_plot('X1(w)', '|X1(w)|', 'fw')
pl.plot(k2, X1)
Out[427]:
In [428]:
pl.figure(figsize=(11,2))
setup_plot('x2(t)', '|x2(t)|', 'time (s)')
pl.plot(t[time_range], x2[time_range]) #snippet between 3s and 3.005s
pl.figure(figsize=(11,2))
setup_plot('X2(w)', '|X2(w)|', 'fw')
pl.plot(k2, X)
Out[428]:
In [429]:
pl.figure(figsize=(11,2))
setup_plot('x(t)', '|x(t)|', 'time (s)')
pl.plot(t[time_range], x[time_range]) #snippet between 3s and 3.005s
pl.figure(figsize=(11,2))
setup_plot('X(w)', '|X(w)|', 'fw')
pl.plot(k2, X)
Out[429]:
In [430]:
IPython.display.Audio(x1, rate=fs)
Out[430]:
In [431]:
IPython.display.Audio(x2, rate=fs)
Out[431]:
In [432]:
IPython.display.Audio(x, rate=fs)
Out[432]:
In [433]:
r = 1/np.sqrt(2)
a1 = [1] + [0]*10
b1 = [1] + [0]*9 + [-1]
In [434]:
print(a1)
print(b1)
In [435]:
z1, p1, k1 = signal.tf2zpk(b1, a1)
p1= []
b1, a1 = signal.zpk2tf(z1, p1, k1)
zplane(z1, p1)
In [436]:
w1, H1 = signal.freqz(b1, a1, whole=True)
H =np.abs(H1)
setup_plot('H(w)', '|H(w)|', 'w')
pl.plot(w1, H)
Out[436]:
In [437]:
y = signal.lfilter(b1, a1, x) #applying the designed filter
IPython.display.Audio(y1, rate=fs)
Out[437]:
Question 2: Comb filter with resonant poles¶
In [438]:
a2 = [1] + [0]*9 + [-0.95]
b2 = [1] + [0]*9 + [-1]
z2, p2, k2 = signal.tf2zpk(b2, a2)
zplane(z2, p2)
print(a2)
print(b2)
In [439]:
w2, H2 = signal.freqz(b2, a2, whole=True)
H2 = np.abs(H2)
setup_plot('H(w)', '|H(w)|', 'w')
pl.plot(w2, H2)
Out[439]:
In [440]:
y2 = signal.lfilter(b2, a2, x) #applying the designed filter
IPython.display.Audio(y2, rate=fs)
Out[440]:
Conclusion¶
Both Comb filters succees in filtering out the noise signal. The All zero comb filter does however weaken the quality of the song, as the guitar in the background becomes quite soft and a bit distorted - thus the all zero filter, succeeds in filtering the noise, but weakens the quality of the music.
The comb filter with resonant poles, also succeeds in filtering out the noise, but in the beginning of the snippet, you can still here some of the noise, until it is filtered out completely once again. Applying this filter, does not effect the quality of the music that much.
In [ ]: